CREATED USING THE RSC LaTeX PCCP ARTICLE TEMPLATE - SEE www.rsc.org/electronicfiles FOR DETAILS 

ARTICLE TYPE www.rsc.org/xxxxxx | XXXXXXXX 



On representing chemical environments 



Albert P. Bartok*", Risi Kondor* and Gabor Csanyi^ 



Received Xth XXXXXXXXXX 20XX, Accepted Xth XXXXXXXXX 20XX 
First published on the web Xth XXXXXXXXXX 200X 
DOI: 10.1039/bOOOOOOx 

T-H We review some recently proposed methods to represent atomic neighbourhood environments and analyse their relative merits 
in terms of their faithfulness and suitability for fitting potential energy surfaces (PES). The crucial properties that such repre- 
sentations (commonly called descriptors) must have are continuity and invariance to the basic symmetries of physics: rotation, 
"q reflection, translation, and permutation of identical atoms. We demonstrate that schemes that initially look quite different are 
^\ specific cases of a general approach, in which a finite set of basis functions with increasing angular wave numbers are used to 
expand the atomic neighbourhood density function. We quantitatively show using the example system of small clusters that this 
expansion needs to be carried to higher and higher wave numbers as the number of neighbours increases in order to obtain a 
faithful representation, and that variants of the descriptors converge at very different rates. 
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1 Introduction 

The appropriate representation of atomic environments in 
terms of a descriptor, a tuple of typically real valued functions 
of the atomic positions, is a crucial ingredient of algorithms 
used in modern computational chemistry and condensed mat- 
ter physics. For example, in structure search applications', 
each configuration depends numerically on the precise ini- 
tial conditions and the path of the search, so it is important 
to be able to identify equivalent structures and detect sim- 
ilarities. In molecular dynamics simulations of phase tran- 
sitions'2, one needs good order parameters that are capable 
of detecting changes in the local order around given atoms. 
When constructing interatomic potentials and fitting potential 
energy surfaces (PES)'^'^, the driving application behind the 
present work, the functional forms depend on components of 
a carefully chosen representation of atomic neighbourhoods, 
e.g. bond lengths, bond angles, etc. "In silico" drug discov- 
eryEIS- and other areas of chemical informatics also rely on 
characterising molecules using similar descriptors (also called 
fingerprints). 

While specifying the position of each atom in a Carte- 
sian coordinate system provides a simple and unequivocal de- 
scription of atomic configurations, it is not directly suitable 
for making comparisons between structures: the list of co- 
ordinates is ordered arbitrarily and two structures might be 
mapped to each other by a rotation, reflection or translation so 
that two different lists of atomic coordinates can, in fact, repre- 
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sent the same or very similar structures. A good representation 
is invariant with respect to permutational, rotational, reflec- 
tional and translational symmetries, while retaining the faith- 
fulness of the Cartesian representation. In particular, a system 
of invariant descriptors , g'2 , . . . , is said to be complete if 
it uniquely determines the atomic environment, up to transla- 
tion, rotation, and reflection. It is said to be over-complete if it 
contains spurious descriptors in the sense that a proper subset 
of {^1,^2, ■ ■ ■ ,?m} is, by itself, complete. If a representation 
is complete, a one-to-one mapping (a bijection) is obtained be- 
tween the genuinely different atomic environments and the set 
of invariants comprising the representation. An over-complete 
representation assigns potentially many distinct descriptors to 
a given atomic structure, but guarantees that genuinely dif- 
ferent atomic structures will never have identical descriptors 
associated with them: the function relating representations to 
atomic structures is a surjection. 

Fitting potential energy surfaces (PESs) of small molecules 
to data generated by first principles electronic structure cal- 
culations has been a mainstay of computational chemistry for 
decadesl^^l^Il jn these applications, the number of atoms which 
comprise the configuration space is relatively small, and, cru- 
cially, fixed. Typically, when modelling the PES of a small 
cluster, pairwise distances - or, equivalently, reciprocal dis- 
tances - are usedR Such descriptions work only for a fixed 
number of atoms in a fixed order, thus they are limited to de- 
scribing a given set of atoms. Even in this case, a seemingly 
new configuration is obtained by just permuting the order of 
atoms, i.e. crucial symmetries may be missing in this frame- 
work. Braams and Bowman*^ remedied this shortcoming by 
using polynomials of basis functions of the pairwise distances, 
selected so that they are invariant to the permutation of identi- 
cal atoms. Computer code is available that generates the per- 
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mutationally invariant polynomials automatically'^ (up to ten 
atoms), but this approach does not allow for varying number 
of atoms in the database of configurations. 

In order to generate interatomic potentials for solids or large 
clusters capable of describing a wide variety of conditions, the 
number of neighbours that contribute to the energetics of an 
atom has to be allowed to vary, with the descriptor remaining 
continuous and differentiable. None of the traditional repre- 
sentations fulfil this criterion. Recently, however, a number 
of new, promising descriptors have been proposed together 
with potential energy surfaces based on them^* ^ ^''^ and at 
this point it is not clear which method of representing atomic 
neighborhoods will prove to be optimal in the long term. In 
order to disentangle this issue from the rather complex details 
of generating first principles data and fitting PES, this paper 
focuses solely on the problem of descriptors. 

The most well known invariants describing atomic neigh- 
bourhoods are the bond-order parameters originally proposed 
by Steinhardt et al. These have been successfully used as or- 
der parameters in studies of nucleation^^, phase transitions^^ 
and glasses.!^ Here we show that the bond-order parameters 
form a subset of a more general set of invariants called the bis- 
pectrum^K The formally infinite array of bispectrum compo- 
nents provides an over-complete system of invariants, and by 
truncating it one obtains representations whose sensitivity can 
be refined at will. We relate the b ispectrum to the representa- 
tion proposed by Behler et a/. 121221^ and show that, together with 
another descriptor set described below, their angular parts are 
all simple polynomials of the same canonical invariants. 

2 Descriptors 

Among the applications mentioned in the introduction, some 
require representing the geometry of an entire molecule, while 
for others one needs to describe the neighbourhood of an atom 
perhaps within a finite cutoff distance. While these two cases 
are closely related, the descriptors for one are not directly suit- 
able for the other, although often the same idea can be used 
to derive representations for either case. Below, we focus on 
representing the neighbour environment of a single atom, but 
for some cases briefly mention easy generalisations that yield 
global molecular descriptors. 

For atomic position vectors {r,! , r,2, . . . , r,7v}, taken rela- 
tive to a central atom /, the symmetric matrix 



r,i ■ r,i r,i • r,2 
ra ■ Til r,'2 • r;2 



Til • TiW 



(1) 



is, according to Weyl,'^an over-complete table of basic invari- 
ants to rotation, reflection and translation. However, E is not a 



suitable descriptor, because permutations of atoms change the 
order of rows and columns. For example, swapping atoms 1 
and 2 results in the transformed matrix 



r/i-r,-2 ra-r,! 
TiAT • r,-2 Tin ■ r/i 



fa ■ r,w 



(2) 



To compare two structures using their Weyl matrices E and 
E', we define a reference distance metric 



-*ref 



minllE-PE'P^I 

p M 



(3) 



where P is a permutation matrix and the minimisation is over 
all possible permutations. This metric is not differentiable at 
locations where the permutation that minimises Q changes. 
It would also be intractable to calculate exactly for large num- 
bers of atoms, but nevertheless we will use this metric to assess 
the faithfulness of other representations. 

One way to generate permutationally invariant continuous 
functions of the Weyl matrix is to compute its eigenvalues, 
indeed, such a descriptor was recently used to fit the atomisa- 
tion energies of a large set of molecules'^. However, the list 
of eigenvalues is very far from being complete, since there are 
only A' eigenvalues, whereas the dimensionality of the config- 
uration space of neighbours is 3A^ — 3 (after the rotational 
symmetries are removed). It is also unclear how to make the 
descriptors based on the eigenspectrum continuous and differ- 
entiable as the number of neighbours varies. Other, continuous 
invariants shown later in this paper are very closely related to 
the elements of E. 

Another straightforward way to compare structures is based 
on pairing the atoms from each and finding the optimal rota- 
tion and reflection that brings the two structures into as close 
an alignment as possible. For each pair of structures {r,;}^^[ 
and {r'yl'^^p it is possible to order the atoms according to 



their distance from the central atom ; 
- and compute 



or use other heuristics 



A(^) 



N 

i 



where R represents an arbitrary rotation and define the dis- 
tance between two configurations as 



: minA(/?). 

R 



(4) 



This distance clearly has all the necessary invariances and 
completeness properties, but, like t/^ef^ it is not suitable for 
parametrising potential energy surfaces, because it is again not 
differentiable: the reordering procedure and the minimisation 
over rotations and reflections introduce cusps. 
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In the field of molecular informatics, one popular descrip- 
tor is based on the histogram of pairwise atomic distances^, 
similar to Valle's crystal fingerprint^^. We will not consider 
it here, because it is unsuitable for fitting PES, as it is clearly 
not complete: e.g. from 6 unordered distance values it is not 
necessarily possible to construct a unique tetrahedron, even 
though the number of degrees of freedom is also 6. 

2.1 Bond-order parameters 

As a first step in deriving a continuous invariant representation 
of atomic environments, we define the local atomic density 
function associated with atom / as 



(5) 



where the index j runs over the neighbours of atom / and r,y 
is the vector from atom / to atom j. In writing (|5]l we have as- 
sumed that only one atomic species is present, but generalising 
to multiple species is straightforward, e.g. by introducing dif- 
ferent weights in the sum, corresponding to different species. 
Determining which neighbours to include in the summation 
can be done by using a simple binary valued, or a smooth 
real valued cutoff function of the interatomic distance, or via 
a more sophisticated procedure e.g. Voronoi analysis. The 
local atomic density is already invariant to permuting neigh- 
bours, because changing the order of the atoms in the neigh- 
bour list only affects the order of the summation. In order to 
simpHfy the following derivation, for now we omit the infor- 
mation on the radial distance to the neighbours, but will show 
later how the radial information can be included. The local 
atomic density function can be expanded in terms of spherical 
harmonics (dropping the atomic index / for clarity): 



l=Qm=-l 



(6) 



where f is the point on the sphere corresponding to the direc- 
tion of the vector r, thus p(f) is the projection of p(r) onto 
the unit sphere. 

The properties of functions defined on the unit sphere are 
related to the group theory of SO(3), the group of three dimen- 
sional rotations. Spherical harmonics form an orthonormal ba- 
sis set for L2{S^), the class of square integrable functions on 
the sphere: 

{Ylm\Yl'm') = ^ll'^mm'^ 

where the inner product of functions / and g is defined as 

{f\8)= /r(fk(f)d£2(f), 

where the surface element dn(f ) can be expressed in terms of 
the polar angles 9 and as 

dn(f) =sin0d0d0. 



and the coefficients c/„, are given by 

Cim ^ {p\Yl„,) =I^i'/m(f,;)- 



(7) 



The quantities Qi,„ introduced by Steinhardt et a/.'^are pro- 
portional to the coefficients c/„,. Averaging the coefficients for 
atom ; provides the atomic order parameters 



' j 



(8) 



where A^, is the number of neighbours of atom /. Furthermore, 
averaging over atoms in the entire system gives a set of global 
order parameters 

where Ni, is the total number of atom pairs included in the 
summation. Both sets are invariant to permutations of atoms 
and translations, but still depend on the orientation of the ref- 
erence frame. However, rotationally invariant combinations 
can be constructed as 



Qi = 



An 

21 + 



m=—l 



1/2 



and 



L 

m\ ,m2:«i3=— / 



I I 



I 



nil ni2 niT, 



(9) 



(10) 



for atoms and 

Qi = 
Wi = 



An 

21 + 



1 1/2 



7 E QlnQn 

1 m=-l 



I 

L 



I I I 

mi m2 niT, 



Qlm\ Qlm2 Qlm2 



for the entke structure. The factor in parentheses is the 
Wigner-3 jm symbol,S9which is zero unless mi +m2+m^ = 0. 

The numbers Q) and Wj are called second-order and third- 
order bond-order parameters, respectively. It is possible to 
normalise W/ such that it does not depend strongly on the 
number of neighbours: 



-1 -3/2 



Im 



=-/ 



For symmetry reasons, only coefficients with I > A have 
non-zero values in environments with cubic symmetry and 
/ > 6 for environments with icosahedral symmetry. Different 
values correspond to crystalline materials with different sym- 
metry, while the global order parameters vanish in disordered 
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phases, such as liquids. Bond-order parameters were origi- 
nally introduced for studying order in liquids and glasses 
but were soon adopted for a wide range of applications. They 
have been used to study the free energy of clusters* ^^ * ^^ ! melt- 
ing of quantum solids '*^, nucleation'^"', as well as to serve as 
reaction coordinates in simulations of phase transitionsl^Sl and 
also to generate interatomic potentials^'. 

Using some basic concepts from representation theory, we 
now prove that the second-order bond-order parameters are in- 
deed rotationally invariant, then we show a more general set 
of third order invariants^', of which the Qs and the Ws are a 
subset. An arbitrary rotation R operating on a spherical har- 
monic function Yi,„ transforms it into a linear combination of 
spherical harmonics with the same / index: 

/ 

m' = -l 

where the D' (R) matrices are known as the Wigner matrices, 
which form the irreducible representations of the three dimen- 
sional rotation group, SO(3). The elements of the Wigner ma- 
trices are given by 

D',„„,iR)^{Yi,„\R\Yi,„,). (11) 
It follows that the rotation operator R acts on the function p as 

t=Qm=-l l=Om=-l 

l=Om=-lm'=-l 
I 

= E E '^'linYlm'^ 

l=Om'=-l 

thus the column vector C/ of expansion coefficients transforms 
under rotation as 

C/^D'(^)c,. 

Making use of the fact that rotations are unitary operations 
on L2{S^), it is possible to show that the matrices D' are uni- 
tary, 

D'V=I, 
leading us to a set of rotational invariants, 

pi=c]ci. (12) 

which we call the rotational power spectrum due to the anal- 
ogy with the familiar power spectrum of ordinary Fourier anal- 
ysis. The power spectrum remains the same under rotations: 



(c|D'^) (d 



We also note that the elements of C/ transform under reflection 
about the origin as 

C/^(-l)'c,, (13) 

thus the power spectrum is invariant to this symmetry oper- 
ation. A comparison with equations (j?]), ([8]) and (j9|l shows 
that the second-order bond-order parameters are related to the 
power spectrum via the simple scaling 



Qi 



4;r 

2Z + 1^' 



The power spectrum is clearly not a complete descriptor for 
a general function /(f) on the sphere, for example consider 
the two different functions 



and 



/l =J22+i'2-2+i33+}'3-3 



/2=F21+i'2-l+i32+F3-2, 



which both have the same power spectrum, p2 =2 and 773 — 
2. However, for the restricted class of functions which are 
sums of delta functions (such as the local atomic density p in 
equation (|5])) with the radial information removed, numerical 
experiments suggest that the power spectrum might be over- 
complete. 

2.2 The bispectrum 

We generalise the concept of the power spectrum to obtain a 
larger set of invarian ts via the coupling of the different angular 
momentum channels'^"''-. Let us consider the direct product 
C/j (E) C/j , which transforms under a rotation as 



C/. (8)C/, ^> 



C/. 



It follows from the representation theory of compact groups 
that the direct product of two irreducible representations of a 
compact group can be decomposed into a direct sum of irre- 
ducible representations of the same groupEl In case of the 
SO(3) group, the direct product of two Wigner matrices can 
be decomposed into a direct sum of Wigner matrices. 



Lm. 



where CmjJi.nn denote the Clebsch-Gordan coefficients or, us- 
ing more compact notation. 



:C C,. 



/=l'l-'2l 



(14) 
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where C'' '^ denote the matrices formed of the Clebsch- 
Gordan coefficients. These are themselves unitary, so the vec- 
tor C''''- (c/j (EiC/j) transforms as 



h+h 

e D' 

l=\h-h\ 



Qh-h 



C/, ®Cu 



(15) 



We now define gi^.i^j by 
'1+/2 

l=\l\-h\ 



i.e. g/, ./2,/ is that part of the RHS in equation ( 15 1 which trans- 
forms under rotation as 

Analogously to the power spectrum, we can now define the 
bispectrum as the collection of scalars 

biih.i =c|g/i./2,/, 
which are invariant to rotations: 



(c,D')'d' 



It follows from equation ( 13 1 that those elements of the bis- 



pectrum where Zi + Z2 + ^ is odd, change sign under reflection 
about the origin. If invariance to reflection is required, we take 
the absolute value of these components or omit them. 
Rewriting the bispectrum formula as 



h 



m=—lmi=—li m2=— /t 



the similarity to the third-order bond-order parameters be- 
comes apparent. Indeed, the Wigner 3 jm-symbols are related 
to the Clebsch-Gordan coefficients through 



h h h 
m\ m2 mj, 



(17) 



and by substituting the spherical harmonics identity F/,,, = 
(-l)"'y,l^ in equation (Isll it foflows that 



(18) 



GL = (-i)'"(ejj* 



Substituting the identities ^TT) and ( fTSj l into the definition ( 10 1 
we obtain 



1 



V2r- 



E (-i)-'"c; 



mi ,m2,m3 =— / 



GUeU(-i)"'(eu)*' 



thus the third-order parameters W/ are seen to be proportional 
to the diagonal elements of the bispectrum, bi j j. Noting that 
Foo = 1, the coefficient coo is the number of neighbours A^, and 
using c'^f^Q^,^ — Si i^dmjni^ the bispectrum elements h =0,1 = 
h are identical to the previously introduced power spectrum 
components: 

m=— /m2=— / 



m=—l 



therefore. 



The first few terms of the power spectrum and bispectrum 
for atom / with three neighbours are shown below, where 0,,^ 
is the angle between the bonds ij and ik and the sums are over 
the neighbours. 



Po 
Pi 

P2 

P3 

PA 



9 

A% 
3 

An 



Ecos% + 3 



— { -V cos^ft,i. + 6 



7 /5 
An 



9 / 35 
\bn \ y 



Ecos4 0,,.,-15Ecos2 0,.y, + 13 



mn^ \ A 



Jk 



"Ecos^e, 



ijk 



jk 



5Ecos0,-^(t 

Jk ) 



hi 



150 / 7 / 5 



13 



7/5 5 

33 U E ^ilk + 7 E ^ilk E COS Qi,k 
^ V jk ^ jk jk 



'iJk 



Jk 



- 4 E cos^ dijk - 2 E cos Oijk + 1 8 

jk jk 



2.3 Radial basis 



Thus far we neglected the distance of neighbouring atoms 
from the central atom by using the unit-sphere projection of 
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the atomic environment. One way to introduce radial informa- 
tion is to complement the spherical harmonics basis in equa- 
tion ^ with radial basis functions 

I 

PW""EE E c„i„,g„ir)Yi„,(r). 

n l=Om=-l 

If the set of radial basis functions is not orthonormal, i.e. 
(gnlgm) = S,„n ^ 5,„„, after obtaining the coefficients c^^,,,, with 

the elements c„/„, are given by 

Cnlm = E ('^ ')n'n '^'n'lm- 
n' 

In practice, when constructing the bispectrum, either c'^^^^^ or 
c„i,„ can be used. 

Rotational invariance must only apply globally, and not to 
each radial basis separately, therefore the angular momentum 
channels corresponding to different radial basis functions need 



to be coupled. So although extending equations ( 12 1 and ( 16 1 
simply as 



Pn.l — E ^nlm^i'lm 

I h h 

E E 



bnJi,l2,l ~ E E E ^nUn^mjn\,m2^nl\''H^nh'"2 



provides a set of invariants describing the three-dimensional 
neighbourhood of the atom, this can easily lead to a poor 
representation. However, if the radial basis functions do not 
sufficiently overlap, the different channels will still not be 
fully coupled, and the representation will have spurious quasi- 
invariances to rotating subsets of atoms at approximately the 
same distance, as illustrated in figure [T] 

To avoid this, it is necessary to choose basis functions that 
are sensitive in a wide range of distances, although this may 
reduce the sensitivity of each radial basis channel, if the func- 
tions are varying very slowly. The fine-tuning of the basis set 
is rather arbitrary, and there is no guarantee that a choice exists 
that is optimal or even satisfactory for all systems of interest. 

We suggest constructing radial functions from cubic and 
higher order polynomials, 0a (r) = (rem — '")"^^/A^a for a = 
1,2,..., «maxj normalised on the range (0, rc ) using 



Na 



(^cut - r)2(«+2)dr 



„2a+5 



' cut 



2a + 5 



The orthonormahsed construction 

gn{r) = ^ W„a0a(r) 
a=l 



(19) 



1 

2 



r2 



?"cut 



Fig. 1 Example of very localised radial basis functions, Gaussians 
in this case. Atoms 1 and 2 at distance r\ and r2 from the centre 
become decoupled as their contribution to the power spectrum or 
bispectrum elements is weighed down by the product gn{,>'\)gn(r2), 
which is rather small for all n. 



guarantees that each radial function returns smoothly to zero 
at the cutoff with both the first and the second derivative being 
zero. The linear combination coefficients are obtained from 
the overlap matrix as 



c r\ ^ ^w, r u \/(5 + 2a)(5 + 2j3) 



5 + a+j8 



W = S-i/2. 




^cut 



Fig. 2 Example of radial basis functions gnir), as defined in 
equation i\9\ for n = 1,2,3,4. 
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Another way to avoid radial decoupling is to define the rota- 
tional invariants in such a way that they couple different radial 
channels explicitly, for example, as 

/ 

Pni.n2J ~ 52 ''inhn'-tnlm 

m=—l 

bni,n2Ji.l2J ~ ^ ^ ^ ''ni/mQ«,mf,m2''«2'l™l^"2'2™2' 

Here, each invariant has contributions from two different ra- 
dial basis channels, and so we ensure that they cannot become 
decoupled, but at the price of increasing the number of invari- 
ants quadratically or even cubically in the number of radial 
basis functions. 

2.4 4-dimen$ional power spectrum and bispectrum 

We now present an alternative to the power spectrum and bis- 
pectrum that does not need the explicit introduction of a ra- 
dial basis set, but still represents atomic neighbourhoods in 
three-dimensional space. We start with projecting the atomic 
neighbourhood density within a cutoff radius rcut onto the sur- 
face of the four-dimensional sphere with radius kq. The 
surface of is defined as the set of points s G M^, where 
S] + ^2 + *3 + *4 = '■q' while the polar angles 0, 9 and 0o of s 
are defined so that 

si = rocos 00 

52 = ro sin Oq cos 

53 = ro sin Oq sin 9 cos (j) 

54 — ro sin Oq sin sin . 

We choose to use the projection 



(j) = arctan(3'/x) 
= arccos(z/|r| 
00 = ^\r\/ro 



where ro > rcut is a parameter, thus rotations in three- 
dimensional space correspond to a subset of rotations in four- 
dimensional space. This projection is somewhat similar to a 
Riemann projection, except in a Riemann projection 0o would 
be defined as 

00 = 2arctan(|r|/2ro) 




implying 



0o« |r|/rofor |r| < ro. 



In contrast to the Riemann projection, our choice of 0o allows 
more sensitive representation of the entire radial range. The 
Umit ro = rcut projects each atom at the cutoff distance to the 
south pole of the 4-dimensional sphere, thus losing all angular 



information. Too large an ro would project all positions onto 
a small surface area of the sphere, requiring a large number of 
basis functions to represent the atomic environment. In prac- 
tice, a large range of ro values works well, in particular, we 
used ro = '"cut- 

To illustrate the procedure. Figure |3] shows the Riemann 
projections for one and two dimensions, which can be easily 
drawn. 




Fig. 3 Two- and three-dimensional Riemann constmctlons that map 
a flat space onto the surface of a sphere. 



An arbitrary function p defined on the surface of a 4D 
sphere can be numerically represente d usin g the hyperspheri- 
cal harmonic functions U^,^{(l),9, 00)'^^'^ 



E E 

j=Oui.m'=—j 



mm mm 



(21) 



which, in fact, correspond to individual matrix components 
of the Wigner (i.e. rotational) matrices, as defined in equa- 
tion (111. In this case the arguments represent a rotation by 0o 



around the vector pointing in the (0,0) direction, which can 
be transformed to the conventional Euler-angles, and j takes 
half-integer values. 

The hyperspherical harmonics form an orthonormal basis, 
thus the expansion coefficients c' , can be calculated via 

^ mm 

cj , ^{U\ |p), 

m m ^ m m i~ ' ' 

where (.|.) denotes the inner product defined on the four- 
dimensional hypersphere: 

{f\g)= / d0osin2 0o / d0sin0 / d(/) 

Jo Jo Jo 

r(0o,0, 0)^(00, 0,0). 

Although the coefficients , have two indices besides /, for 
each /■ it is convenient to collect them into a single vector c-'. 
Similarly to the three-dimensional case, rotations act on the 
hyperspherical harmonic functions as 



m J m ] m J m 1 m2 f«2 ^2 "^2 

m'^m2 
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where the matrix elements , , are given by 



m^m\m2m2 



{U-', \R\U', ) 



Hence the rotation R acting on p transforms the coefficient 
vectors c-' according to 



The unitary R' matrices are the SO(4) analogues of the 
Wigner-matrices D' of the SO(3) case above, and it can be 
shown that the direct product of the four-dimensional rotation 
matrices decomposes according to 



r.7l,;2^ 



.71 +72 
i=l7l-./2 



H 



71 J2 



which is the 4-dimensional analogue of equation (14i. The 
coupling constants H^' -'^, or Clebsch-Gordan coefficients of 
SO(4) areEam 

j[mim'^j2fn2ni2 m^mi,m2 m',m'^.m2' 

The rest of the derivation continues analogously to the 3D 
case, and finally we arrive at the expression for the SO(4) bis- 
pectrum elements 

71 72 7 . . ... 

■^7lJ2,7= 52 52 ^m'!mi%2^m'lm\,m'2^ 

m'j ,mi =— m2,m2=— 72 in',m=—j 

and the 4D power spectrum can be constructed as 



c\ c', . 

mm I mm 



The 4D bispectrum is invariant to rotations of four- 
dimensional space, which include three-dimensional rotations. 
However, there are additional rotations, associated with the 
third polar angle Qq, which, in our case, represents the radial 
information. In order to eliminate the unphysical invariance 
with respect to rotations along the third polar angle, we mod- 
ify the atomic density as 

p,(r) = 5(0)+^5(r-r„-), 

7 

i.e. we add the central atom as a fixed reference point. The 
resulting invariants Bj^j^.j have only three indices, but con- 
tain both radial and angular information, and have the required 
symmetry properties. There are no adjustable parameters in 



the definition of the invariants, apart from the radial cutoff, 
Tcut, and the projection parameter tq. 

The number of components in the truncated representation 
depends on the band limit jmax in the expansion ( 21 1. For sym- 
metry reasons, the bispectrum components with non-integer 
ji + ji + J change sign under reflection and, because of this 
reason, we omitted them. Just as in the 3D case, the represen- 
tation is probably over-complete, i.e. most of the bispectrum 
components are redundant. To reduce the number of redun- 
dant elements, we only used the 'diagonal' components, i.e. 
ji — j2. Table[T]shows the number of bispectrum elements for 
increasing band limit values. 

Table 1 Number of compontents in the full and diagonal bispectrum 
as a function of the band limit jmax- 
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2.5 Parrinello-Behler descriptor 

We include in the tests below the descriptor suggested by Par- 
rinello and Behler^ using the parameters published recentlyE3 
(and henceforth termed PB). The two- and three-body symme- 
try functions (in their terminology) are, 

gL = E = ^'^P [-na{nj-Rsaf] fcinj) 



and 



Gta^2'-^" L (l+A„cos0,,-,)?°exp[-T7„(r2 +4 + r2 

'X.fc{rij)fc{rik)fc{rjk), 



where the cutoff function is defined as 

1 



/.('•) = 



cos 







/2 for r < rcut 
for r > rcut 



Different values of the parameters ri,Rs,^,X can be used to 
generate an arbitrary number of invariants. 

2.6 Angular Fourier Series 

Notice that the angular part of the power spectrum, bispec- 



trum (section 2.2 1 and the descriptors defined by Parrinello 



and Behler (section 2.5 1 are all simple polynomials of the 
canonical set Y.jk^'^^"' (^ijk for integer to, which in turn are 
sums of powers of the basic invariants of Weyl. We include 
in the tests in the next section a further descriptor set, which 
we call the Angular Fourier Series (AFS) descriptor, formed 
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Fig. 4 Radial basis functions in tiie Parrinello-Beiiler (PB) type 
descriptors.'^ 

by the orthogonal polynomials of the basic invariants, which 
are conveniently chosen as the Chebyshev-polynomials Ti{x), 
as 

r/(cos0) =cos(/0), 

and incorporate the radial information using the basis func- 
tions defined in equation ( [T9] l, leading to 

^^n] = Y^8n{rij)gn{rik)co^{lOiik)- 
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Fig. 5 Examples of the angular basis functions for /max = 4 of the 
AFS descriptor. 



3 Reconstruction experiments 

Recall that the elements r,j • r,^ of E defined in equation ([TJi 
are an over-complete set of basic invariants, which, in the case 
of atoms scattered on the surface of a unit sphere (|r,j| = 1), 
are the cosines of the bond angles, Oij^. Thus, the angular part 
of the above descriptors are permutationally invariant func- 
tions of the basic invariants in E, and depending on the actual 
number of descriptor elements used, they may form a com- 
plete or over-complete representation of the atomic environ- 
ment. There are (A^^ + A^) /2 basic invariants, thus we need at 
most this many independent descriptor elements, whereas the 
number of independent degrees of freedom in the neighbour- 
hood configuration is 3A^ — 3, so it remains unclear how many 
of the descriptor elements are actually needed in order to re- 
construct an atomic environment up to an arbitrary rotation 
and permutation of the neighbour atoms. Obtaining exact 
results regarding completeness would be difficult, especially 
since the algebraically dependency relationships between the 
descriptor elements is unknown. 

It is worth emphasising that severely incomplete descrip- 
tors will lead to poor potential energy surfaces, because en- 
tire manifolds of atomic configurations with widely varying 
true energies will be represented by the same descriptor and 
hence assigned the same energy in the fitted PES, resulting 
in many degenerate modes. As the mappings are highly non- 
linear, in general we could not find an analytical route to invert 
the descriptor transformations and thus to analytically inves- 
tigate the completeness properties. However, it is possible to 
conduct numerical experiments in which we compare the de- 
scriptors of a fixed target with that of a candidate structure and 
minimise the difference with respect to the atomic coordinates 
of the candidate. The purpose of these experiments is to de- 
termine if a representation is complete or not, and in the latter 
case to characterise the degree of its faithfulness. 

The global minimum of the descriptor difference between 
the target and the candidate is the zero value and is always 
a manifold due to the symmetries built into the descriptors, 
but for an incomplete descriptor, many inequivalent structures 
will also appear equivalent, thus enlarging the dimensionality 
of the global minimum manifold. Further, it can be expected 
that the descriptor difference function has a number of local 
minima. 

In our experiments we tried to recover a given target struc- 
ture after randomising its atomic coordinates. For each n 
(4 < n < 19) we used 10 different Si„ clusters, obtained from 
a tight-binding^** molecular dynamics trajectory at a temper- 
ature of 2000 K. For each cluster, we selected an atom as the 
origin, randomised its neighbours within the radial cutoff of 
the descriptors by some amount, and then tried to reconstruct 
the original structure by minimising the magnitude of the dif- 
ference between the descriptors of the fixed target and the can- 
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didate as the atomic positions were varied. We repeated this 
procedure five times for each target configuration and then se- 
lected the reconstruction which ended up closest to the target. 
In some cases, such as large n and large randomisation, it was 
possible for our procedure to find only a local minimum dur- 
ing the minimisation, these were easily identified (by noting 
the small gradient of the objective function while the value of 
the objective function was not small) and discarded. 

In order to assess the success of the reconstruction proce- 
dure we employed the reference measures defined in equa- 
tions ([3]) and However, in some cases it was difficult or 
impossible to find the right rotation R in whereas d^^f in 
([3]l proved reliable. For each d^ef, an initial P was generated 
by ordering the atoms according to their distances from the 
central atom, then the optimal permutation was found using a 
simple random search in the space of permutations. 

We minimised the difference between the target and can- 
didate descriptors in the space of atomic coordinates using 
the Conjugate Gradients algorithm, stopping the minimisa- 
tion if either the gradient or the reference distance dygf be- 
came smaller than 10^^ and 10^^, respectively. In order to 
ensure that structures deemed non-equivalent by t/^f > 10^^ 
were genuinely different, we cross-checked them by noting 
the value of A from equation (j4]) and also employing the 
atomic fingerprints suggested by Oganov and Valle ^^^. To give 
a sense of the typical magnitude of d^ef measure, the actual 
difference in terms of atomic distances between two example 
structures is shown in Figure |6j 




Fig. 6 Two Sig clusters that differ by rfj-ef = 4. 1 A^. Atoms ; and /' at 
the origin, i.e. the centre of rotations, is coloured black. In terms of 
Parinello-Behler type descriptors, the difference La(Gia — Gfa)^ 
between the atomic environments of i (left) and i' (right) is 6 • 10^^. 
The bond lengths are shown, in angstroms. 



band limit of the angular descriptors (corresponding to our 
Wx or jmax) was ^max = 16 and only the values ^ = 1,2,4, 16 
are used. 

Figure |7] shows the quality of reconstruction based on the 
PB, AFS, 3D power spectrum, 3D bispectrum and 4D bispec- 
trum as given by the average reference distance d^^f achieved. 
The general trend is the same for all descriptors: as the num- 
ber of neighbours increases, the average dt^f increases, and 
thus the faithfullness of the recontruction decreases. Noting 
that the stopping criteron for the reconstruction process was 
c/ref < 10^^, larger randomisation of the initial atomic coor- 
dinates (bottom panel) reveals the poor representation power 
for all descriptors for n > 10, and the neighbour configuration 
becomes impossible to determine from the descriptor. 

The poor quality of representation is partly attributable to 
the decrease in sensitivity to the positions of atoms near the 
cutoff. For example, Figure[8]shows two Sig clusters for which 
the all the descriptors failed to lead to perfect reconstructions 
(resulting in the observed peak on Figure]?]). The atom marked 
A in the figure is within the 6 A cutoff, but close to it. In or- 
der to separate out this effect, we repeated the reconstruction 
experiments with a radial cutoff of 9 A (omitting the PB de- 
scriptor now since there is no published optimal parameter set 
for this cutoff). The results are shown in Figure|9]for the larger 
initial randomisation. The peak near n = 8 is now absent, and 
the transition from faithful reconstruction (for n < 9 for the 
3D power spectrum and AFS, and for n < 12 for the 4D bis- 
pectrum) to failure for larger n is much clearer. 

Since all the descriptors are likely to be over-complete when 
the infinite series of the basis set expansion is not truncated, 
the reconstruction quality is expected to increase with increas- 



ing descriptor length. To verify this. Figure 10 shows the re 



construction quality of the AFS descriptor for varying trun- 
cations of the angular part of the basis set. The representa- 
tion becomes monotonically better for higher angular resolu- 
tions. However, this comes at the price of introducing ever 
more highly oscillating basis functions, which might be less 
and less suitable for fitting generally smooth potential energy 
surfaces. 



In the first set of recontruction experiments, in order to pro- 
vide a fair comparison, the truncation of the formally infinite 
set of descriptors were chosen in such a way that the finite de- 
scriptors had roughly equal number of components: 5 1 in total 
for the bispectrum and PB descriptors and 50 for the AFS and 
power spectrum. This corresponds to a truncation of the 4-D 
bispectrum with jmax = 5, the PB descriptor with its published 
parameters and the AFS and power spectrum using Z^ax = 9 
and Hjnax = 5. We note that in case of the PB descriptor the 



To verify that the above results are not affected by arte- 
facts of the minimisation procedure, e.g. getting stuck. Fig- 
ure 1 1 shows the convergence of the reference measure d^et 



during a minimisation as well as the corresponding conver- 
gence of the target function (the difference between the target 
and candidate descriptors). There was no difficulty in converg- 
ing the target function to zero (the global minimum) for any 
of the descriptors, while the reference similarity converged to 
a nonzero value for incomplete descriptors. 
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4 Conclusion 

In this paper we discussed a number of approaches to represent 
atomic environments such that the representation is invariant 
to global rotations, reflections, and permutations of atoms. We 
showed that the bond-order parameters introduced by Stein- 
hardt et al. are equivalent to certain elements of the SO(3) 
power spectrum and bispectrum. To incorporate radial infor- 
mation and therefore provide a full description of the atomic 
environment we reviewed the construction of the SO(4) power 
spectrum and bispectrum as an alternative to introducing ex- 
plicit radial basis functions. We also demonstrated that these 
constructs, as well as the descriptors suggested by Parrinello 
and Behler, in fact, use very similar terms and form part of 
a general family based on a Fourier series expansion of bond 
angles. In practice, when the Fourier series is truncated, the 
faithfulness of the descriptors decreases as the number of 
neighbours increases. The descriptors discussed in this pa- 
per are thus suitable for comparing atomic environments, with 
a tuneable tradeoff between the size (and thus computational 
cost) of evaluating the descriptor and its faithfulness in terms 
of its ability to represent the atomic environment uniquely up 
to symmetries. With typically used parameters, however, all 
descriptors failed for Si clusters with more than 13 atoms. The 
relative performance of the various representations in the con- 
text of fitting potential energy surfaces remains to be tested. 
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Fig. 7 Average difference after reconstructing structures 
randomised by 0.2 A (top) and 1 .6 A (bottom) as a function of the 
number of atoms in the cluster. The cutoff was 6 A. Different hnes 
correspond to different descriptors: Parrinello-Behler (PS), Angular 
Fourier Series (AFS), bispectrum (BS) and power spectrum (PS). 
The two versions of the bispectrum differ in the handling of the 
radial degrees of freedom. 
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Fig. 9 Average difference after reconstructing structures 
randomised by 1 .6 A as a function of the number of atoms in the 
cluster. The cutoff was 9 A. The line types corresponding to the 
different descriptors are the same as in Figure |7] 
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Fig. 10 Average difference after reconstructing structures 
randomised by 1.6 A as a function of the number of atoms in the 
cluster using the AFS descriptor with a radial cutoff of 9 A. The 

different curves correspond to different number of components of 
the AFS descriptor, achieved by varying the truncation of the 
angular expansion. 
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Fig. 11 Convergence of reference distance measure during the 
reconstruction procedure. The inset shows the value of the 
minimisation target approaching zero, i.e. descriptor equivalence. 
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